Prostorni podaci su namenjeni skladištenju, manipulaciji i analizi informacija o nekim prostornim fenomenima. Na primer, izmerena temperatura u jednoj tački može da predstavlja jedan prostorni podatak. Prostorni podatak je najčešće reprezentacija jednog prostornog fenomena u toj tački. U ovom slučaju to je temperatura, koja je prostorno-vremenski promenljiva. Prema tome, prostorni podaci opisuju prostorne fenomene i kao takvi, prilagođeni su prirodi i karakteristikama prostornih fenomena:
sf paketRazvoj novih tehnologija i softvera otvorenog koda za rad sa prostornim podacima, kao i standarda u ovoj oblasti, zahteva jednostavniji i efikasniji pristup nad podacima. R paket sp je razvijen da proširi R funkcionalnosti klasama i metodama namenjenim prostornim podacima. Postoje opravdani razlozi za prelazak na „najsavremeniju“ sf klasu, koja omogućava interoperabilnost i efikasnije rukovanje prostornim podacima. Paket sf je R paket za čitanje, čuvanje, rukovanje i manipulaiciju prostornim podacima u R okruženju, putem pristupa koji se naziva simple feature. Simple feature je međunarodni standard za predstavljanje i kodiranje prostornih podataka, predstavljenih putem geometrije tipa tačka, linija i poligon (ISO, 2004). Najvažniji razlog za upotrebu ovog pristupa je bolja manipulacija geometrijom prostornih podataka i jednostavan način za čuvanje geometrije u podacima, tj. jednom redu data.frame-a , čime se u potpunosti mogu iskoristiti funkcionalnosti familije paketa tidyverse.
Tekst zadatka
Vrednosti emisije zagađenja po polutantima SO2, NOx, PM10, PM2.5, NMVOC i NH3 prostorizovati po izvorima zagađenja (tačkastih, linijskih i površinskih kategorija). Vrednosti sumirati po ćelijama grida za područje teritorije Republike Srbije.
ESRI shapefile.geopackage formatu za razmenu prostornih podataka, kao i gridovane podatke po kategorijama kao posebne lejere.Podaci su dobijeni za sledeće kategorije:
Učitavanje tačkastih izvora zagađivača i transformacija:
## Reading layer `1A1a_ep' from data source `D:\R_projects\Nauka_R\Slides\R\Spatial data vezbe\Data\1A1a_ep.shp' using driver `ESRI Shapefile'
## Simple feature collection with 21 features and 14 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 353634.2 ymin: 4795824 xmax: 587823.5 ymax: 5105985
## projected CRS: WGS 84 / UTM zone 34N
## Simple feature collection with 21 features and 14 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 353634.2 ymin: 4795824 xmax: 587823.5 ymax: 5105985
## projected CRS: WGS 84 / UTM zone 34N
## First 10 features:
## NOx SO2 PM10 PM2_5 NMVOC NH3 Unit Quantty Type
## 1 11.5790 0.6673455 0.09684276 0.03978916 0.08637947 0 NA NA NA
## 2 72.5170 4.1750000 1.06799662 0.43880094 0.54097761 0 NA NA NA
## 3 284.5313 36.2251000 3.57381777 1.46835164 2.12260659 0 NA NA NA
## 4 27.4530 1.1820151 0.19795276 0.08133158 0.20479968 0 NA NA NA
## 5 112.1674 27.8750000 2.25393515 0.92605991 0.83677002 0 NA NA NA
## 6 21.2126 1.2580000 0.42127192 0.17308530 0.15824623 0 NA NA NA
## 7 23.1560 0.9121350 0.21193984 0.08707837 0.17274401 0 NA NA NA
## 8 39.2910 96.0010000 18.70001630 7.68315603 0.29311129 0 NA NA NA
## 9 51.9012 0.3656000 0.39069000 0.16052030 0.38718352 0 NA NA NA
## 10 23.6490 2.5205763 0.17792041 0.07310102 0.17642180 0 NA NA NA
## X___10 Region Area
## 1 Beogradske elektrane, TO Banovo brdo Beogradski region Grad Beograd
## 2 Beogradske elektrane, TO Cerak Beogradski region Grad Beograd
## 3 Beogradske elektrane, TO Novi Beograd Beogradski region Grad Beograd
## 4 Beogradske elektrane, TO Miljakovac Beogradski region Grad Beograd
## 5 Beogradske elektrane, TO Dunav Beogradski region Grad Beograd
## 6 Beogradske elektrane, TO Medaković Beogradski region Grad Beograd
## 7 Beogradske elektrane, TO Voždovac Beogradski region Grad Beograd
## 8 Beogradske elektrane, TO Zemun Beogradski region Grad Beograd
## 9 Beogradske elektrane, TO Konjarnik Beogradski region Grad Beograd
## 10 Beogradske elektrane, TO Mirijevo Beogradski region Grad Beograd
## Mncplty Locatin geometry
## 1 Beograd-Čukarica Beograd (Čukarica) POINT (453732.9 4958584)
## 2 Beograd-Čukarica Beograd (Čukarica) POINT (453736.5 4953799)
## 3 Beograd-Novi Beograd Beograd (Novi Beograd) POINT (454145.1 4961030)
## 4 Beograd-Rakovica Beograd (Rakovica) POINT (456831.2 4955267)
## 5 Beograd-Stari Grad Beograd (Stari Grad) POINT (456269 4963981)
## 6 Beograd-Voždovac Beograd (Voždovac) POINT (459853.8 4958226)
## 7 Beograd-Voždovac Beograd (Voždovac) POINT (459553.2 4955717)
## 8 Beograd-Zemun Beograd (Zemun) POINT (451897.5 4965798)
## 9 Beograd-Zvezdara Beograd (Zvezdara) POINT (453732.9 4958584)
## 10 Beograd-Zvezdara Beograd (Zvezdara) POINT (463136.6 4959595)
Učitavanje linijskih izvora zagađivača i transformacija:
## Reading layer `1A3bi_R' from data source `D:\R_projects\Nauka_R\Slides\R\Spatial data vezbe\Data\1A3bi_R.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1440 features and 23 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 333615.1 ymin: 4656734 xmax: 657521.3 ymax: 5114098
## projected CRS: WGS 84 / UTM zone 34N
## Simple feature collection with 1440 features and 23 fields
## geometry type: LINESTRING
## dimension: XY
## bbox: xmin: 333615.1 ymin: 4656734 xmax: 657521.3 ymax: 5114098
## projected CRS: WGS 84 / UTM zone 34N
## First 10 features:
## ID is_PGDS PGDS_2015_ Name Kategrj Broj_pt Oznak_d Smer Oznk_pc Oznk_zv
## 1 1 0 1504.233 IIA IIA 100 10001 O 10001 1301
## 2 2 1 1859.667 IIA IIA 100 10002 O 1302 10002
## 3 3 0 4787.406 IIA IIA 100 10003 O 10002 102
## 4 4 1 4787.371 IIA IIA 100 10004o1 O 102 102.1
## 5 5 0 6670.213 IIA IIA 100 10004o2 D 102.1 10003
## 6 6 0 6669.752 IIA IIA 100 10004o3 L 10003 102.1
## 7 7 0 6670.213 IIA IIA 100 10005o1 O 10003 10003.1
## 8 8 0 6666.552 IIA IIA 100 10005o2 D 10003.1 10003.2
## 9 9 0 6670.213 IIA IIA 100 10005o3 L 10003.2 10003.1
## 10 10 0 4273.142 IIA IIA 100 10005o4 O 10003.2 1103
## Duzina Stacinz Nazv_pc Nzv_zvr
## 1 4.679 4.679 granica MA<U+00D0>/SRB (Horgoš 2) Horgoš (put A1)
## 2 10.439 15.118 Horgoš (Kanjiža) Backi Vinogradi
## 3 3.5 18.618 Backi Vinogradi petlja Subotica sever
## 4 2.349 20.967 petlja Subotica sever <Null>
## 5 9.265 30.232 <Null> Subotica (centar)
## 6 9.244 <Null> Subotica (centar) <Null>
## 7 0.983 31.215 Subotica (centar) <Null>
## 8 2.394 33.609 <Null> <Null>
## 9 2.402 <Null> <Null> <Null>
## 10 2.934 36.543 <Null> Subotica (B.Topola)
## Komentr PGDS_2015 NOx SO2 PM10 PM2_5 NMVOC
## 1 <Null> 1504.233 0.8676208 0.03817683 0.05566276 0.04222431 0.3049724
## 2 <Null> 1377.000 0.7942347 0.03494771 0.05095463 0.03865285 0.2791768
## 3 <Null> 4787.406 2.7613098 0.12150244 0.17715357 0.13438407 0.9706120
## 4 <Null> 6673.000 3.8488946 0.16935807 0.24692826 0.18731333 1.3529026
## 5 <Null> 6670.213 3.8472872 0.16928734 0.24682514 0.18723510 1.3523376
## 6 <Null> 6669.752 3.8470213 0.16927564 0.24680808 0.18722216 1.3522442
## 7 <Null> 6670.213 3.8472872 0.16928734 0.24682514 0.18723510 1.3523376
## 8 <Null> 6666.552 3.8451755 0.16919442 0.24668966 0.18713233 1.3515953
## 9 <Null> 6670.213 3.8472872 0.16928734 0.24682514 0.18723510 1.3523376
## 10 <Null> 4273.142 2.4646895 0.10845063 0.15812371 0.11994851 0.8663487
## NH3 PGDS geometry
## 1 0.03728672 1504.233 LINESTRING (418315.4 511096...
## 2 0.03413289 1377.000 LINESTRING (410263.6 510721...
## 3 0.11866955 4787.406 LINESTRING (406835.8 510659...
## 4 0.16540940 6673.000 LINESTRING (404627.4 510609...
## 5 0.16534032 6670.213 LINESTRING (398266.3 510462...
## 6 0.16532890 6669.752 LINESTRING (404627.1 510609...
## 7 0.16534032 6670.213 LINESTRING (397392 5104191,...
## 8 0.16524957 6666.552 LINESTRING (397840.2 510196...
## 9 0.16534032 6670.213 LINESTRING (397394.8 510419...
## 10 0.10592205 4273.142 LINESTRING (396709.4 509936...
Učitavanje površinskih izvora zagađivača - pomoćni podaci i transformacija:
## Reading layer `CLC18_RS' from data source `D:\R_projects\Nauka_R\Slides\R\Spatial data vezbe\Data\CLC18_RS.shp' using driver `ESRI Shapefile'
## Simple feature collection with 34794 features and 6 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 331043.4 ymin: 4674468 xmax: 663570 ymax: 5117030
## projected CRS: ETRS89 / UTM zone 34N
## Simple feature collection with 34794 features and 6 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 331043.4 ymin: 4674468 xmax: 663570 ymax: 5117030
## projected CRS: ETRS89 / UTM zone 34N
## First 10 features:
## ID CODE_18 Area_Ha Remark SHAPE_Leng SHAPE_Area
## 1 RS_1 111 57.65009 <NA> 6456.549 576500.9
## 2 RS_2 111 37.62057 <NA> 5540.471 376205.7
## 3 RS_3 111 35.45603 <NA> 3725.036 354560.3
## 4 RS_4 111 25.14260 <NA> 4178.823 251426.0
## 5 RS_5 112 42.43190 <NA> 4099.389 424319.0
## 6 RS_6 112 67.68868 <NA> 4377.618 676886.8
## 7 RS_7 112 28.49615 <NA> 2736.777 284961.5
## 8 RS_8 112 25.46822 <NA> 2175.572 254682.2
## 9 RS_9 112 39.81418 <NA> 3961.963 398141.8
## 10 RS_10 112 34.70737 <NA> 2637.301 347073.7
## geometry
## 1 POLYGON ((574615.8 4796652,...
## 2 POLYGON ((572794.7 4796959,...
## 3 POLYGON ((457109.9 4962932,...
## 4 POLYGON ((397083.5 5106354,...
## 5 POLYGON ((555111.5 4676318,...
## 6 POLYGON ((555062.8 4679080,...
## 7 POLYGON ((560340.8 4680704,...
## 8 POLYGON ((553657.7 4681341,...
## 9 POLYGON ((559394.7 4682410,...
## 10 POLYGON ((553051 4682493, 5...
Učitavanje grida i transformacija u koordinatni sistem u projekciji (UTM projekcija, 34N zona / WGS84 elipsoid):
## Reading layer `Grid_Serbia_005' from data source `D:\R_projects\Nauka_R\Slides\R\Spatial data vezbe\Data\Grid_Serbia_0.05deg.gpkg' using driver `GPKG'
## Simple feature collection with 3769 features and 6 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 18.8 ymin: 42.2 xmax: 23 ymax: 46.2
## geographic CRS: WGS 84
## Simple feature collection with 3769 features and 6 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 18.8 ymin: 42.2 xmax: 23 ymax: 46.2
## geographic CRS: WGS 84
## First 10 features:
## left top right bottom id value geom
## 1 18.80 45.95 18.85 45.90 6 1 POLYGON ((18.8 45.95, 18.85...
## 2 18.80 45.90 18.85 45.85 7 1 POLYGON ((18.8 45.9, 18.85 ...
## 3 18.80 45.85 18.85 45.80 8 1 POLYGON ((18.8 45.85, 18.85...
## 4 18.85 45.95 18.90 45.90 96 1 POLYGON ((18.85 45.95, 18.9...
## 5 18.85 45.90 18.90 45.85 97 1 POLYGON ((18.85 45.9, 18.9 ...
## 6 18.85 45.85 18.90 45.80 98 1 POLYGON ((18.85 45.85, 18.9...
## 7 18.85 45.80 18.90 45.75 99 1 POLYGON ((18.85 45.8, 18.9 ...
## 8 18.85 45.75 18.90 45.70 100 1 POLYGON ((18.85 45.75, 18.9...
## 9 18.85 45.65 18.90 45.60 102 1 POLYGON ((18.85 45.65, 18.9...
## 10 18.85 45.60 18.90 45.55 103 1 POLYGON ((18.85 45.6, 18.9 ...
vars <- c("NOx","SO2", "PM10", "PM2_5", "NMVOC", "NH3")
sf.1A1a_ep %<>% dplyr::select(vars)
names(sf.1A1a_ep)## [1] "NOx" "SO2" "PM10" "PM2_5" "NMVOC" "NH3" "geometry"
Prostorizacija - sumiranje vrednosti tačkastih zagađivača po ćelijama grida:
p.1A1a_ep <- gridWGS84 %>%
sf::st_join(sf.1A1a_ep) %>%
dplyr::group_by(id) %>%
dplyr::summarize(NOx = sum(NOx, na.rm = TRUE),
SO2 = sum(SO2, na.rm = TRUE),
PM10 = sum(PM10, na.rm = TRUE),
PM2_5 = sum(PM2_5, na.rm = TRUE),
NMVOC = sum(NMVOC, na.rm = TRUE),
NH3 = sum(NH3, na.rm = TRUE)) %>%
dplyr::mutate(id = as.numeric(id))Kontrola:
sum.sf.1A1a_ep <- sf.1A1a_ep %>%
sf::st_drop_geometry() %>%
dplyr::select(., vars) %>%
apply(., 2, sum) %>%
t(.) %>%
as.data.frame() %>%
dplyr::mutate_if(is.numeric, round, 2)
sum.p.1A1a_ep <- p.1A1a_ep %>%
sf::st_drop_geometry() %>%
dplyr::select(., vars) %>%
apply(., 2, sum) %>%
t(.) %>%
as.data.frame() %>%
dplyr::mutate_if(is.numeric, round, 2)
data.frame(sum = c("prostorizovano", "ukupno", "razlika"), rbind(sum.p.1A1a_ep, sum.sf.1A1a_ep, data.frame(sum.p.1A1a_ep - sum.sf.1A1a_ep))) %>%
datatable(caption = "Tabela 1: Sumarne vrednosti razlika nakon prostorizacije") Vizuelizacija:
map.1A1a_ep <- p.1A1a_ep
map.1A1a_ep$Spatialised <- 0
map.1A1a_ep$Spatialised[map.1A1a_ep$NOx != 0 | map.1A1a_ep$SO2 != 0 | map.1A1a_ep$PM10 != 0 | map.1A1a_ep$PM2_5 != 0 | map.1A1a_ep$NMVOC != 0 | map.1A1a_ep$NH3 != 0] <- 1 # Zadatak: resiti sa mutate i case_when
mapview(map.1A1a_ep, zcol = "Spatialised", layer.name = "Spatialised 1A1a_ep")Mreža puteva po kategorijama puteva:
ggplot()+
geom_sf(data = sf.1A3bi_R, aes(colour = Name))+
scale_colour_manual(name = "Kategorija: ", values = c("red", "orange"))+
labs(title = "Mreža puteva",
subtitle = "Klasifikacija po kategorijama puteva",
caption = "Republika Srbija")+
theme_bw()+
theme(legend.position = "bottom")Mreža puteva po vrednosti protoka saobraćaja:
ggplot()+
geom_sf(data = sf.1A3bi_R, aes(size = PGDS_2015/10000), colour = "black")+
scale_size_identity()+
labs(title = "Mreža puteva",
subtitle = "Klasifikacija po vrednosti protoka saobraćaja",
caption = "Republika Srbija")+
theme_bw()+
theme(legend.position = "bottom")Sumiranje ukupnih vrednosti po polutantima:
sum.sf.1A3bi_R <- sf.1A3bi_R %>%
sf::st_drop_geometry() %>%
dplyr::select(., vars) %>%
apply(., 2, sum) %>%
t(.) %>%
as.data.frame() %>%
dplyr::mutate_if(is.numeric, round, 2)
sum.sf.1A3bi_R %>% datatable()Presek sa gridom:
sf.1A3bi_R %<>% dplyr::mutate(Length = st_length(.),
PGDSL = PGDS_2015 * Length)
sf.1A3bi_R.int <- st_intersection(sf.1A3bi_R, gridWGS84) %>%
dplyr::select(., PGDS_2015, PGDSL, Length) %>%
mutate(Length_int = st_length(.),
PGDS_int = (PGDS_2015/Length)*Length_int) %>%
dplyr::select(., Length, PGDS_int)
sum_PGDS <- sum(sf.1A3bi_R.int$PGDS_int)
sf.1A3bi_R.int <- sf.1A3bi_R.int %>%
dplyr::mutate(NOx = ((sum.sf.1A3bi_R$NOx/sum_PGDS)*PGDS_int),
SO2 = ((sum.sf.1A3bi_R$SO2/sum_PGDS)*PGDS_int),
PM10 = ((sum.sf.1A3bi_R$PM10/sum_PGDS)*PGDS_int),
PM2_5 = ((sum.sf.1A3bi_R$PM2_5/sum_PGDS)*PGDS_int),
NMVOC = ((sum.sf.1A3bi_R$NMVOC/sum_PGDS)*PGDS_int),
NH3 = ((sum.sf.1A3bi_R$NH3/sum_PGDS)*PGDS_int)) %>%
dplyr::select(vars)
# Kontrola
sum.sf.1A3bi_R.int <- sf.1A3bi_R.int %>%
sf::st_drop_geometry() %>%
dplyr::select(., vars) %>%
apply(., 2, sum) %>%
t(.) %>%
as.data.frame() %>%
dplyr::mutate_if(is.numeric, round, 2)
sum.sf.1A3bi_R## NOx SO2 PM10 PM2_5 NMVOC NH3
## 1 3099.65 136.39 198.86 150.85 1089.54 133.21
## NOx SO2 PM10 PM2_5 NMVOC NH3
## 1 3099.65 136.39 198.86 150.85 1089.54 133.21
Prostorizacija - sumiranje vrednosti linijskih zagađivača po ćelijama grida:
p.1A3bi_R <- gridWGS84 %>%
st_join(sf.1A3bi_R.int, join = st_contains) %>%
mutate(id = as.numeric(id)) %>%
group_by(id) %>%
summarize(NOx = sum(NOx, na.rm = TRUE),
SO2 = sum(SO2, na.rm = TRUE),
PM10 = sum(PM10, na.rm = TRUE),
PM2_5 = sum(PM2_5, na.rm = TRUE),
NMVOC = sum(NMVOC, na.rm = TRUE),
NH3 = sum(NH3, na.rm = TRUE)) %>%
mutate(id = as.numeric(id))
p.1A3bi_R[,vars] <- drop_units(p.1A3bi_R[,vars])Kontrola:
sum.sf.1A3bi_R <- sf.1A3bi_R %>%
sf::st_drop_geometry() %>%
dplyr::select(., vars) %>%
apply(., 2, sum) %>%
t(.) %>%
as.data.frame() %>%
dplyr::mutate_if(is.numeric, round, 2)
sum.p.1A3bi_R <- p.1A3bi_R %>%
sf::st_drop_geometry() %>%
dplyr::select(., vars) %>%
apply(., 2, sum) %>%
t(.) %>%
as.data.frame() %>%
dplyr::mutate_if(is.numeric, round, 2)
data.frame(sum = c("prostorizovano", "ukupno", "razlika"), rbind(sum.p.1A3bi_R, sum.sf.1A3bi_R, data.frame(sum.p.1A3bi_R - sum.sf.1A3bi_R))) %>%
datatable(caption = "Tabela 2: Sumarne vrednosti razlika nakon prostorizacije") Vizuelizacija:
map.1A3bi_R <- p.1A3bi_R
map.1A3bi_R$Spatialised <- 0
#map.1A3bi_R[,vars] <- units::drop_units(map.1A3bi_R[,vars])
map.1A3bi_R$Spatialised[map.1A3bi_R$NOx != 0 | map.1A3bi_R$SO2 != 0 | map.1A3bi_R$PM10 != 0 | map.1A3bi_R$PM2_5 != 0 | map.1A3bi_R$NMVOC != 0 | map.1A3bi_R$NH3 != 0] <- 1
mapview(map.1A3bi_R, zcol = "Spatialised", layer.name = "Spatialised 1A3bi_R")Corine Land Cover karta (zemljišni pokrivač):
ggplot()+
geom_sf(data = CLC_12_18, aes(fill = CODE_18))+
#scale_colour_manual(name = "Kategorija: ", values = c("red", "orange"))+
labs(title = "Corine Land Cover karta",
subtitle = "Klasifikacija po kategorijama zemljišnog pokrivača",
caption = "Republika Srbija")+
theme_bw()+
theme(legend.position = "right")Ukupne vrednosti po polutantima:
sum.sf.3Da1 <- data.frame(NOx = 8.33, SO2 = 0, PM10 = 0, PM2_5 = 0, NMVOC = 0, NH3 = 18028.55)
sum.sf.3Da1## NOx SO2 PM10 PM2_5 NMVOC NH3
## 1 8.33 0 0 0 0 18028.55
Izdvajanje poljoprivrednog zemljišta (CLC classes) i presek sa gridom:
sf_clc18_polj <- subset(CLC_12_18, CODE_18 == "211" | CODE_18 == "221" | CODE_18 == "222" | CODE_18 == "242" | CODE_18 == "243") # CLC agricultural areas
sf_clc18_polj[,vars] <- NAPresek sa gridom:
sf_clc18_polj.int <- st_intersection(sf_clc18_polj, gridWGS84) %>%
dplyr::mutate(Area = st_area(.)) # racunanje povrsina poligona
sf_clc18_polj.int %<>% dplyr::select(., vars, Area) Interpolacija vrednosti po polutantima srazmerno površini entiteta:
sum_Area <- sum(sf_clc18_polj.int$Area)
sf.3Da1 <- sf_clc18_polj.int %>%
mutate(NOx = ((sum.sf.3Da1$NOx/sum_Area)*Area),
SO2 = ((sum.sf.3Da1$SO2/sum_Area)*Area),
PM10 = ((sum.sf.3Da1$PM10/sum_Area)*Area),
PM2_5 = ((sum.sf.3Da1$PM2_5/sum_Area)*Area),
NMVOC = ((sum.sf.3Da1$NMVOC/sum_Area)*Area),
NH3 = ((sum.sf.3Da1$NH3/sum_Area)*Area))
sf.3Da1 %<>% dplyr::select(vars)Prostorizacija - sumiranje vrednosti površinskih zagađivača po ćelijama grida:
p.3Da1 <- gridWGS84 %>%
st_join(sf.3Da1, join = st_contains) %>%
group_by(id) %>%
summarize(NOx = sum(NOx, na.rm = TRUE),
SO2 = sum(SO2, na.rm = TRUE),
PM10 = sum(PM10, na.rm = TRUE),
PM2_5 = sum(PM2_5, na.rm = TRUE),
NMVOC = sum(NMVOC, na.rm = TRUE),
NH3 = sum(NH3, na.rm = TRUE)) %>%
mutate(id = as.numeric(id))
p.3Da1[, vars] <- drop_units(p.3Da1[, vars])Kontrola:
sum.p.3Da1 <- p.3Da1 %>%
sf::st_drop_geometry() %>%
dplyr::select(., vars) %>%
apply(., 2, sum) %>%
t(.) %>%
as.data.frame() %>%
dplyr::mutate_if(is.numeric, round, 2)
data.frame(sum = c("prostorizovano", "ukupno", "razlika"), rbind(sum.p.3Da1, sum.sf.3Da1, data.frame(sum.p.3Da1 - sum.sf.3Da1))) %>%
datatable(caption = "Tabela 2: Sumarne vrednosti razlika nakon prostorizacije") Vizuelizacija:
map.3Da1 <- p.3Da1
map.3Da1$Spatialised <- 0
#map.3Da1[,vars] <- units::drop_units(map.3Da1[,vars])
map.3Da1$Spatialised[map.3Da1$NOx != 0 | map.3Da1$SO2 != 0 | map.3Da1$PM10 != 0 | map.3Da1$PM2_5 != 0 | map.1A3bi_R$NMVOC != 0 | map.3Da1$NH3 != 0] <- 1
mapview(map.3Da1, zcol = "Spatialised", layer.name = "Spatialised 1A3bi_R")data.spat.list <- list(p.1A1a_ep, p.1A3bi_R, p.3Da1)
sf_data <- data.spat.list[[1]]
for(i in 2:length(data.spat.list)){
sf_data <- st_join(sf_data, data.spat.list[[i]], join = st_equals) %>%
group_by(id.x) %>%
summarize(NOx = sum(NOx.x, NOx.y),
SO2 = sum(SO2.x, SO2.y),
PM10 = sum(PM10.x, PM10.y),
PM2_5 = sum(PM2_5.x, PM2_5.y),
NMVOC = sum(NMVOC.x, NMVOC.y),
NH3 = sum(NH3.x + NH3.y)) %>%
mutate(id = id.x) %>%
select(id, NOx, SO2, PM10, PM2_5, NMVOC, NH3)
print(paste("NOx:",sum(sf_data$NOx)))
print(paste("SO2:",sum(sf_data$SO2)))
print(paste("PM10:",sum(sf_data$PM10)))
print(paste("PM2_5:",sum(sf_data$PM2_5)))
print(paste("NMVOC:",sum(sf_data$NMVOC)))
print(paste("NH3:",sum(sf_data$NH3)))
}## [1] "NOx: 4339.7786"
## [1] "SO2: 1910.84306396663"
## [1] "PM10: 323.042992264955"
## [1] "PM2_5: 201.872271321156"
## [1] "NMVOC: 1098.79137281681"
## [1] "NH3: 133.21"
## [1] "NOx: 4348.1086"
## [1] "SO2: 1910.84306396663"
## [1] "PM10: 323.042992264955"
## [1] "PM2_5: 201.872271321156"
## [1] "NMVOC: 1098.79137281681"
## [1] "NH3: 18161.76"
## Simple feature collection with 3769 features and 7 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: 18.8 ymin: 42.2 xmax: 23 ymax: 46.2
## geographic CRS: WGS 84
## # A tibble: 3,769 x 8
## id NOx SO2 PM10 PM2_5 NMVOC NH3 geom
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <POLYGON [arc_degree]>
## 1 6 0. 0 0 0 0 0 ((18.8 45.95, 18.85 45.9~
## 2 7 4.25e-5 0 0 0 0 0.0920 ((18.8 45.9, 18.85 45.9,~
## 3 8 1.84e-4 0 0 0 0 0.398 ((18.8 45.85, 18.85 45.8~
## 4 96 1.23e-4 0 0 0 0 0.266 ((18.85 45.95, 18.9 45.9~
## 5 97 2.16e-1 0.00941 0.0137 0.0104 0.0751 4.31 ((18.85 45.9, 18.9 45.9,~
## 6 98 1.08e-1 0.00470 0.00686 0.00520 0.0376 1.70 ((18.85 45.85, 18.9 45.8~
## 7 99 0. 0 0 0 0 0 ((18.85 45.8, 18.9 45.8,~
## 8 100 0. 0 0 0 0 0 ((18.85 45.75, 18.9 45.7~
## 9 102 0. 0 0 0 0 0 ((18.85 45.65, 18.9 45.6~
## 10 103 0. 0 0 0 0 0 ((18.85 45.6, 18.9 45.6,~
## # ... with 3,759 more rows
classes.NOx <- classIntervals(sf_data$NOx, n = 12, style = "fisher")
classes.SO2 <- classIntervals(sf_data$SO2, n = 12, style = "fisher")
classes.PM10 <- classIntervals(sf_data$PM10, n = 12, style = "fisher")
classes.PM2_5 <- classIntervals(sf_data$PM2_5, n = 12, style = "fisher")
classes.NMVOC <- classIntervals(sf_data$NMVOC, n = 12, style = "fisher")
classes.NH3 <- classIntervals(sf_data$NH3, n = 12, style = "fisher")
sf_data <- sf_data %>%
mutate(percent_class_NOx = cut(NOx, classes.NOx$brks, include.lowest = T,dig.lab=5),
percent_class_SO2 = cut(SO2, classes.SO2$brks, include.lowest = T,dig.lab=7),
percent_class_PM10 = cut(PM10, classes.PM10$brks, include.lowest = T,dig.lab=5),
percent_class_PM2_5 = cut(PM2_5, classes.PM2_5$brks, include.lowest = T,dig.lab=5),
percent_class_NMVOC = cut(NMVOC, classes.NMVOC$brks, include.lowest = T,dig.lab=5),
percent_class_NH3 = cut(NH3, classes.NH3$brks, include.lowest = T,dig.lab=5)
)
pal1 <- viridisLite::viridis(12, direction = -1)a<-ggplot() +
geom_sf(data = sf_data,
aes(fill = percent_class_NOx)) +
scale_fill_manual(values = pal1,
name = "NOx [t]") +
labs(x = NULL, y = NULL,
title = "Pollutant inventory spatialization - NOx",
subtitle = "Spatial resolution 0.05°x0.05°",
caption = "UBFCE (2020)") +
theme(line = element_blank(),
#axis.text = element_blank(),
axis.title = element_blank(),
#legend.position = "None", ###################### legend
panel.background = element_blank()) +
coord_sf(datum = sf::st_crs(4326))b<-ggplot() +
geom_sf(data = sf_data,
aes(fill = percent_class_SO2)) +
scale_fill_manual(values = pal1,
name = "SO2 [t]") +
labs(x = NULL, y = NULL,
title = "Pollutant inventory spatialization - SO2",
subtitle = "Spatial resolution 0.05°x0.05°",
caption = "UBFCE (2020)") +
theme(line = element_blank(),
#axis.text = element_blank(),
axis.title = element_blank(),
#legend.position = "None", ###################### legend
panel.background = element_blank()) +
coord_sf(datum = sf::st_crs(4326))c<-ggplot() +
geom_sf(data = sf_data,
aes(fill = percent_class_PM10)) +
scale_fill_manual(values = pal1,
name = "PM10 [t]") +
labs(x = NULL, y = NULL,
title = "Pollutant inventory spatialization - PM10",
subtitle = "Spatial resolution 0.05°x0.05°",
caption = "UBFCE (2020)") +
theme(line = element_blank(),
#axis.text = element_blank(),
axis.title = element_blank(),
#legend.position = "None", ###################### legend
panel.background = element_blank()) +
coord_sf(datum = sf::st_crs(4326))d<-ggplot() +
geom_sf(data = sf_data,
aes(fill = percent_class_PM2_5)) +
scale_fill_manual(values = pal1,
name = "PM2.5 [t]") +
labs(x = NULL, y = NULL,
title = "Pollutant inventory spatialization - PM2.5",
subtitle = "Spatial resolution 0.05°x0.05°",
caption = "UBFCE (2020)") +
theme(line = element_blank(),
#axis.text = element_blank(),
axis.title = element_blank(),
#legend.position = "None", ###################### legend
panel.background = element_blank()) +
coord_sf(datum = sf::st_crs(4326))e<-ggplot() +
geom_sf(data = sf_data,
aes(fill = percent_class_NMVOC)) +
scale_fill_manual(values = pal1,
name = "NMVOC [t]") +
labs(x = NULL, y = NULL,
title = "Pollutant inventory spatialization - NMVOC",
subtitle = "Spatial resolution 0.05°x0.05°",
caption = "UBFCE (2020)") +
theme(line = element_blank(),
#axis.text = element_blank(),
axis.title = element_blank(),
#legend.position = "None", ###################### legend
panel.background = element_blank()) +
coord_sf(datum = sf::st_crs(4326))f<-ggplot() +
geom_sf(data = sf_data,
aes(fill = percent_class_NH3)) +
scale_fill_manual(values = pal1,
name = "NH3 [t]") +
labs(x = NULL, y = NULL,
title = "Pollutant inventory spatialization - NH3",
subtitle = "Spatial resolution 0.05°x0.05°",
caption = "UBFCE (2020)") +
theme(line = element_blank(),
#axis.text = element_blank(),
axis.title = element_blank(),
#legend.position = "None", ###################### legend
panel.background = element_blank()) +
coord_sf(datum = sf::st_crs(4326))